#include <qua/fixed.h>
#include <qua/cast.h>
#include <qua/text.h>
#include <qua/error.h>
#include <cstdio>
#include <cmath>

namespace qua
{
    const uint64_t FIXED_POW10[] = {
            1uLL,                   // 10^0
            10uLL,                  // 10^1
            100uLL,                 // 10^2
            1000uLL,                // 10^3
            10000uLL,               // 10^4
            100000uLL,              // 10^5
            1000000uLL,             // 10^6
            10000000uLL,            // 10^7
            100000000uLL,           // 10^8
            1000000000uLL,          // 10^9
            10000000000uLL,         // 10^10
            100000000000uLL,        // 10^11
            1000000000000uLL,       // 10^12
            10000000000000uLL,      // 10^13
            100000000000000uLL,     // 10^14
            1000000000000000uLL,    // 10^15
            10000000000000000uLL,   // 10^16
            100000000000000000uLL,  // 10^17
            1000000000000000000uLL, // 10^18
    };

    const uint64_t FIXED_HALF_POW10[] = {
            0,                      // 10^0  / 2
            5uLL,                   // 10^1  / 2
            50uLL,                  // 10^2  / 2
            500uLL,                 // 10^3  / 2
            5000uLL,                // 10^4  / 2
            50000uLL,               // 10^5  / 2
            500000uLL,              // 10^6  / 2
            5000000uLL,             // 10^7  / 2
            50000000uLL,            // 10^8  / 2
            500000000uLL,           // 10^9  / 2
            5000000000uLL,          // 10^10 / 2
            50000000000uLL,         // 10^11 / 2
            500000000000uLL,        // 10^12 / 2
            5000000000000uLL,       // 10^13 / 2
            50000000000000uLL,      // 10^14 / 2
            500000000000000uLL,     // 10^15 / 2
            5000000000000000uLL,    // 10^16 / 2
            50000000000000000uLL,   // 10^17 / 2
            500000000000000000uLL,  // 10^18 / 2
    };

    const uint64_t FIXED_DIVIDER = 1000000000000000000uLL;

    const uint64_t FIXED_HALFDIV = FIXED_DIVIDER / 2uLL;

    const uint64_t FIXED_SQRTDIV = 1000000000uLL;

    const uint64_t FIXED_EPSILON  = 1000uLL;

    const uint64_t FIXED_EPSFLOAT = 1000000000000uLL;

    const int FIXED_MAX_PRECISION = 18;

	fixed::fixed()
	    : m_int(),
	      m_frac(),
	      m_prec()
	{
	}

	void fixed::normalize()
	{
        static const uint64_t delta = FIXED_EPSILON >> 1;
	    uint64_t rem = m_frac % FIXED_EPSILON;
	    if( rem )
	    {
            if( rem > delta )
                m_frac += FIXED_EPSILON - rem;
            else
                m_frac -= rem;
	    }
	}

	void fixed::normfloat()
	{
	    static const uint64_t delta = FIXED_EPSFLOAT >> 1;
	    uint64_t rem = m_frac % FIXED_EPSFLOAT;
	    if( rem )
	    {
	        if( rem > delta )
	            m_frac += FIXED_EPSFLOAT - rem;
	        else
	            m_frac -= rem;
	    }
	}

    void fixed::set_prec( int prec )
    {
        if( prec > FIXED_MAX_PRECISION )
            prec = FIXED_MAX_PRECISION;
        m_prec = prec;
    }

	fixed::fixed( double arg, int prec )
	    : m_int( cast<int64_t>( std::floor( arg ) ) ),
	      m_frac( cast<uint64_t>( std::floor( (arg - m_int) * FIXED_DIVIDER + 0.5 ) ) ),
	      m_prec( prec )
	{
	    normalize();
	}

	fixed& fixed::operator = ( double arg )
	{
	    m_int = cast<int64_t>( std::floor( arg ) );
	    m_frac = cast<uint64_t>( arg - m_int );
        normalize();
		return *this;
	}

    fixed::fixed( float arg, int prec )
        : m_int( cast<int64_t>( std::floor( arg ) ) ),
          m_frac( cast<uint64_t>( std::floor( (arg - m_int) * FIXED_DIVIDER + 0.5f ) ) ),
          m_prec( prec )
    {
        normfloat();
    }

    fixed& fixed::operator = ( float arg )
    {
        m_int = static_cast<int64_t>( std::floor( arg ) );
        m_frac = static_cast<uint64_t>( arg - m_int );
        normfloat();
        return *this;
    }

	fixed::fixed( int64_t arg, int prec )
	    : m_int( arg ),
	      m_frac( 0 ),
	      m_prec( prec )
	{
	}

	fixed& fixed::operator = ( int64_t arg )
	{
	    m_int = arg;
	    m_frac = 0;
		return *this;
	}

	fixed::fixed( uint64_t arg, int prec )
	    : m_int( cast<int64_t>( arg ) ),
	      m_frac( 0 ),
	      m_prec( prec )
	{
	}

	fixed& fixed::operator = ( uint64_t arg )
	{
	    m_int = cast<int64_t>( arg );
	    m_frac = 0;
		return *this;
	}

    fixed::fixed( int32_t arg, int prec )
        : m_int( arg ),
          m_frac( 0 ),
          m_prec( prec )
    {
    }

    fixed& fixed::operator = ( int32_t arg )
    {
        m_int = arg;
        m_frac = 0;
        return *this;
    }

    fixed::fixed( uint32_t arg, int prec )
        : m_int( arg ),
          m_frac( 0 ),
          m_prec( prec )
    {
    }

    fixed& fixed::operator = ( uint32_t arg )
    {
        m_int  = arg;
        m_frac = 0;
        return *this;
    }

	fixed::fixed( int16_t arg, int prec )
	    : m_int( arg ),
	      m_frac( 0 ),
	      m_prec( prec )
	{
	}

	fixed& fixed::operator = ( int16_t arg )
	{
	    m_int = arg;
	    m_frac = 0;
		return *this;
	}

	fixed::fixed( uint16_t arg, int prec )
	    : m_int( cast<int64_t>( arg ) ),
	      m_frac( 0 ),
	      m_prec( prec )
	{
	}

	fixed& fixed::operator = ( uint16_t arg )
	{
	    m_int = cast<int64_t>( arg );
	    m_frac = 0;
		return *this;
	}

    fixed::fixed( int8_t arg, int prec )
        : m_int( arg ),
          m_frac( 0 ),
          m_prec( prec )
    {
    }

    fixed& fixed::operator = ( int8_t arg )
    {
        m_int = arg;
        m_frac = 0;
        return *this;
    }

    fixed::fixed( uint8_t arg, int prec )
        : m_int( arg ),
          m_frac( 0 ),
          m_prec( prec )
    {
    }

    fixed& fixed::operator = ( uint8_t arg )
    {
        m_int  = arg;
        m_frac = 0;
        return *this;
    }

    uint64_t fixed::divider()
    {
        return FIXED_DIVIDER;
    }

    double fixed::multiplier()
    {
        static const double fmultiplier = 1.0 / static_cast< double >( FIXED_DIVIDER );
        return fmultiplier;
    }

    int fixed::max_prec()
    {
        return FIXED_MAX_PRECISION;
    }

    fixed::operator double () const
    {
        return m_int + frac();
    }

    fixed::operator int64_t () const
    {
        return m_frac < FIXED_HALFDIV ? m_int : m_int + 1;
    }

    fixed::operator bool () const
    {
        return m_int && m_frac;
    }

    fixed::operator float () const
    {
        return static_cast<float>( operator double () );
    }

    fixed::operator int32_t () const
    {
        return cast<int32_t>( m_int );
    }

    fixed::operator uint32_t () const
    {
        return cast<uint32_t>( m_int );
    }

    fixed::operator int16_t() const
    {
        return cast<int16_t>( m_int );
    }

    fixed::operator uint16_t() const
    {
        return cast<uint16_t>( m_int );
    }

    fixed::operator int8_t() const
    {
        return cast<int8_t>( m_int );
    }

    fixed::operator uint8_t() const
    {
        return cast<uint8_t>( m_int );
    }

	fixed& fixed::operator += ( const fixed& arg )
	{
	    m_int += arg.m_int;
	    m_frac += arg.m_frac;
	    if( m_frac >= FIXED_DIVIDER )
	        ++m_int, m_frac -= FIXED_DIVIDER;
	    m_prec = std::max( m_prec, arg.m_prec );
		return *this;
	}

	fixed& fixed::operator -= ( const fixed& arg )
	{
	    m_int -= arg.m_int;
	    if( m_frac >= arg.m_frac )
	        m_frac -= arg.m_frac;
	    else
	        --m_int, m_frac = FIXED_DIVIDER - (arg.m_frac - m_frac);
        m_prec = std::max( m_prec, arg.m_prec );
		return *this;
	}

	fixed& fixed::operator *= ( const fixed& arg )
	{
		return ( *this = *this * arg );
	}

	fixed& fixed::operator /= ( const fixed& arg )
	{
		return ( *this = *this / arg );
	}

	fixed fixed::operator + ( const fixed& arg ) const
	{
	    fixed res = *this;
		return res += arg;
	}

	fixed fixed::operator - ( const fixed& arg ) const
	{
	    fixed res = *this;
		return res -= arg;
	}

    fixed fixed::operator - () const
    {
        fixed res( -m_int, m_prec );
        if( m_frac )
        {
            --res.m_int;
            res.m_frac = FIXED_DIVIDER - m_frac;
        }
        return res;
    }

	fixed fixed::operator * ( const fixed& arg ) const
	{
	    static const int64_t SIGNBIT = 1LL << 63;
	    static const uint64_t LIM63 = (1uLL << 63) - 1;
	    static const uint64_t LIM18 = LIM63 / FIXED_DIVIDER;
	    static const uint64_t LIM09 = LIM63 / FIXED_SQRTDIV;
	    static const char* ERRMSG = "Multiplication overflow for too big fixed point values.";

	    fixed left = abs();
	    fixed right = arg.abs();
	    uint64_t a = static_cast<uint64_t>( left.m_int );
        uint64_t b = left.m_frac;
	    uint64_t c = static_cast<uint64_t>( right.m_int );
	    uint64_t d = right.m_frac;
	    uint64_t A[4][4], B[4], C[4];
	    B[0] = a / FIXED_SQRTDIV, B[1] = a % FIXED_SQRTDIV, B[2] = b / FIXED_SQRTDIV, B[3] = b % FIXED_SQRTDIV;
	    C[0] = c / FIXED_SQRTDIV, C[1] = c % FIXED_SQRTDIV, C[2] = d / FIXED_SQRTDIV, C[3] = d % FIXED_SQRTDIV;

	    if( C[0] && B[0] > LIM63 / C[0] )
            throw error( ERRMSG, __FILE__, __LINE__ );

	    for( int i = 0; i < 4; ++i )
	    {
	        for( int j = 0; j < 4; ++j )
	        {
	            A[i][j] = B[i] * C[j];
	            if( A[i][j] > LIM63 )
	                throw error( ERRMSG, __FILE__, __LINE__ );
	        }
	    }

	    if( A[0][0] > LIM18 )
	        throw error( ERRMSG, __FILE__, __LINE__ );

	    if( A[0][1] > LIM09 || A[1][0] > LIM09 )
	        throw error( ERRMSG, __FILE__, __LINE__ );

	    uint64_t E = A[0][0] * FIXED_DIVIDER;
        if( E > LIM63 )
            throw error( ERRMSG, __FILE__, __LINE__ );

	    E += ( A[0][1] + A[1][0] ) * FIXED_SQRTDIV;
	    if( E > LIM63 )
	        throw error( ERRMSG, __FILE__, __LINE__ );

	    E += A[0][2] + A[1][1] + A[2][0];
        if( E > LIM63 )
            throw error( ERRMSG, __FILE__, __LINE__ );

        uint64_t D = A[0][3] + A[1][2] + A[2][1] + A[3][0];
        if( D > LIM63 )
            throw error( ERRMSG, __FILE__, __LINE__ );

	    E += D / FIXED_SQRTDIV;
        if( E > LIM63 )
            throw error( ERRMSG, __FILE__, __LINE__ );

        uint64_t F = ( D % FIXED_SQRTDIV ) * FIXED_SQRTDIV + A[1][3] + A[2][2] + A[3][1] + ( A[2][3] + A[3][2] ) / FIXED_SQRTDIV;
        if( F % 10 == 9 ) ++F; // correction of 999... at end
        if( F >= FIXED_DIVIDER )
        {
            E += F / FIXED_DIVIDER;
            F  = F % FIXED_DIVIDER;
        }

        fixed res( static_cast<int64_t>( E ), std::max( m_prec, arg.m_prec ) );
        res.m_frac = F;

        int64_t negative = (m_int & SIGNBIT) ^ (arg.m_int & SIGNBIT);

        return negative ? -res : res;
	}

	fixed fixed::operator / ( const fixed& arg ) const
	{
		return *this * arg.reverse();
	}

    fixed fixed::reverse() const
    {
        fixed res( 0LL, m_prec );
        fixed arg = abs();
        uint64_t a = static_cast<uint64_t>( arg.m_int );
        uint64_t b = arg.m_frac;
        if( !a )
        {
            if( !m_frac )
                throw error( "Division by zero value with fixed point.", __FILE__, __LINE__ );
            res.m_int  = static_cast<int64_t>( FIXED_DIVIDER / b );
            res.m_frac = (FIXED_DIVIDER % b) * FIXED_DIVIDER;
        }
        else
        {
            int n = 1;
            for( ; n <= 18 && a >= FIXED_POW10[n]; ++n );
            if( n > 18 ) res.m_frac = 0;
            else
            {
                uint64_t u = FIXED_DIVIDER;
                uint64_t v = a * FIXED_POW10[ FIXED_MAX_PRECISION - n ] + b / FIXED_POW10[ n ];
                res.m_frac = u / v * FIXED_POW10[ FIXED_MAX_PRECISION - n ];
                for( ++n; n <= 18; ++n )
                {
                    u = u % v * 10;
                    if( !u ) break;
                    res.m_frac += u / v * FIXED_POW10[ FIXED_MAX_PRECISION - n ];
                }
            }
        }
        return m_int < 0 ? -res : res;
    }
}
